VE_SD_SchoolConnectedness
Leigh Gayle Goetschius
Date of Last Edit: 27 October, 2020
Import Data
#load('/Volumes/csmonk/csmonk-lab/Data/Fragile_Families_Longitudinal/Leigh_FFCWS/FFCW_core_ffid_transfer.rdata')
FullData = read.spss('/Volumes/csmonk/csmonk-lab/Data/Fragile_Families_Longitudinal/Leigh_FFCWS/FFCW_core_ffid_transfer.sav',to.data.frame=TRUE, use.value.labels = FALSE)Preparing the Data
The composite scores for violence exposure and social deprivation are calculated exactly how they’re done in Hein et al., 2020 in SCAN.
Creating Some Initial Flags
###Creating a variable to indicate whether child lived with mother most of the time at each time point
FullData$age3cust <- ifelse(FullData$m3a2 == 1 |FullData$m3a2 == 2, "1", "0")
FullData$age5cust <- ifelse(FullData$m4a2 == 1 |FullData$m4a2 == 2, "1", "0")
FullData$age9cust <- ifelse(FullData$m5a2 == 1 |FullData$m5a2 == 2, "1", "0")###Creating a variable to indicate whether mom was currently married, living with, or romantically involved with baby’s father at age 9 (since m5c4 doesn’t seem to be recognized) ####FFCWS documentation indicates that they would assign a 1 / yes for the following responses to a4 - 1(married), 4(cohabiting or living together), 5(romantically involved but living apart), -1(refused), -2(don’t know) ####Based on conversations with Colter Mitchell about refusing itself being of interest, will instead only include1, 4, or 5
FullData$con.m5c4 = ifelse(FullData$m5a4 == 1 | FullData$m5a4 == 4 | FullData$m5a4 == 5, "1", "0")###Creating a variable to indicate whether mom was currently in any relationship at each time point
FullData$age3rel = ifelse(FullData$m3d6 == 1 |FullData$m3e2 == 1, "1", "0")
FullData$age5rel = ifelse(FullData$m4d5 == 1 |FullData$m4e2 == 1, "1", "0")
FullData$age9rel = ifelse(FullData$con.m5c4 == 1 |FullData$m5d2 == 1, "1", "0")###Creating a flag to indicate that the mom said she was not in a relationship - with BF or CP
FullData$age3.norel = ifelse(FullData$m3d6 == 2 & FullData$m3e2 == 2, "1", "0")
FullData$age5.norel = ifelse(FullData$m4d5 == 2 & FullData$m4e2 == 2, "1", "0")
FullData$age9.norel = ifelse(FullData$con.m5c4 == 2 & FullData$m5d2 == 2, "1", "0")###Creating a flag to indicate that the mom refused to answer question about relationship with BF
FullData$age3.refusebf = ifelse(FullData$m3a4 == -1, "1", "0")
FullData$age5.refusebf = ifelse(FullData$m4a4 == -1, "1", "0")
FullData$age9.refusebf = ifelse(FullData$m5a4 == -1, "1", "0")###Creating a flag to indicate that the mom refused to answer question about relationship with a CP
FullData$age3.refusecp = ifelse(FullData$m3e2 == -1, "1", "0")
FullData$age5.refusecp = ifelse(FullData$m4e2 == -1, "1", "0")
FullData$age9.refusecp = ifelse(FullData$m5d2 == -1, "1", "0")Recode for Missing
-9 (not in wave), -6 (skip), -3(missing), -2(don’t know), -1 (refuse)
is.na(FullData) <- FullData == -9
is.na(FullData) <- FullData == -8
is.na(FullData) <- FullData == -6
is.na(FullData) <- FullData == -5
is.na(FullData) <- FullData == -3
is.na(FullData) <- FullData == -2
is.na(FullData) <- FullData == -1Select Variables
Violence Exposure and Social Deprivation Composite Scores The items used in the violence exposure and social deprivation composites come from Hein et al., 2020, SCAN are are the same as those in Goetschius et al., 2020, DCN, Goetschius et al., 2020, JAMA Network Open, and Peckins et al., 2019, PNE. I haven’t repeated where they come from here, but their origins can be found in those papers. Additionally, specific papers that T.Hein selected items based on should be referenced in the code below.
School Connectedness 9 & 15 Items for this come from the Panel Study on Income Dynamics Child Development Supplement (PSID-CDS) (The Panel Study of Income Dynamics Child Development Supplement: User Guide for CDS-III, 2010)
Age 9 Internalizing and Externalizing Symptoms Items for both of these came from the Child Behavior Checklist (Achenbach & Edelbrock, 1983). The internalizing items came from the anxious/depressed, withdrawn/depressed, and somatic symptoms subscales and the externalizing items came from the rule-breaking and aggressive behavior subscales.
Positive Adolescent Function Items from this come from the EPOCH Measure of Adolescent Wellbeing (Kern et al., 2016)
Adolescent Internalizing Symptoms Items for this come from: - Center for Epidemiologic Studies Depression Scale (CES-D) (Radloff, 1977) - Brief Symptom Inventory 18 (BSI-18) (Derogatis & Savitz, 2000)
Adolescent Externalizing Symptoms Items for this come from: - Delinquency: National Longitudinal Study of Adolescent Health (Add Health - (Harris, 2013) - Dickman’s Impulsivity scale (Dickman, 1990) - Binary substance use variables (alcohol-2+ drinks, cigarettes, marijuana, illegal drugs other than marijuana, prescription drugs w/out prescription)
Covariates Most of the covariates are pretty self-explanatory, but one requires a bit more information. To get at income I used the poverty ratio that’s constructed in the Fragile Families dataset. Because we use data from most of the waves and I don’t want to have multiple covariates representing income in my SEM for model degrees of freedom purposes, I calculate the average poverty ratio and use that as the covariate. According to the Fragile Families website, the poverty ratio is the ratio of total household income, as defined in c_1hhinc, to the official poverty thresholds, designated by the U.S. Census Bureau. The thresholds in c_1povca vary by family composition and year. At each wave, we used the poverty thresholds for the year preceding the interview. We calculated separate thresholds based on mother and father reports of household size and composition. However, calculations for married/cohabiting mothers and fathers rely on mother reports of household size and composition. A small number of missing values (don’t know, refused) were treated as 0 in household membership counts.
Primary Caregiver Relationship Used this to determine who the primary caregiver generally was for the reported data at the study waves included in the violence exposure and social deprivation composites as well as the externalizing and internalizing at age 9.
Data.ThreatDepVars = FullData %>%
dplyr::select( m1b2, cm1edu,cm2md_case_con,cm2md_case_lib,m2b17a,m2b17b,m2b17c,m2b17d,m2b17e,m2b17f,IH3j10,IH3j11,IH3j13,IH3j14,IH3j15,IH3j16,IH3j17,IH3j18,IH3j19,IH3j3,IH3j4,IH3j6,IH3j7,IH3j8,IH3j9,IH3k1a,IH3k1b,IH3k1c,IH3k1d,IH3k1e,IH3k2a,IH3k2b,IH3k2c,IH3k2d,IH3k2e,IH3l1,IH3l2,IH3l3,IH3l4,IH3l5,IH3l6,IH3l7,m3a2,m3a4,m3d6,m3d9n1,m3d9p,m3d9p1,m3d10,m3d10d,m3d1f,m3d7a,m3d7b,m3d7c,m3d7d,m3d7e,m3d7f,m3d7g,m3d7h,m3d7i,m3d7j,m3d7k,m3d7l,m3d7n,m3d7n1,m3d7p,m3d7p1,m3e2,m3e21b,m3e21f,m3e23a,m3e23b,m3e23c,m3e23d,m3e23e,m3e23f,m3e23g,m3e23h,m3e23i,m3e23j,m3e23k,m3e23l,m3e24,m3h3,m3h3a,m3h4,m3h5,m3h6,m3i6c, m3i6f,m3i23f,m3j36j,m3j36k,m3k26a, m3k26b,IH3r10,IH3s2,IH3s3,IH3s4,IH5g3,IH5g4,IH5g6,IH5g7,IH5g8,IH5g9,IH5g10,IH5g11,IH5g13,IH5g14,IH5g15,IH5g16,IH5g17,IH5g18,IH5g19,IH5h1,IH5h2,IH5h3,IH5h4,IH5h5,IH5h6,IH5h7,m4a2,m4a4,m4d1b,m4d1f,m4d5,m4d7a,m4d7b,m4d7c,m4d7d,m4d7e,m4d7f,m4d7g,m4d7h,m4d7i,m4d7j,m4d7o,m4d7p,m4d10,m4d10a,m4d10e,m4e2,m4e21b,m4e21f,m4e23a,m4e23b,m4e23c,m4e23d,m4e23e,m4e23f,m4e23g,m4e23h,m4e23i,m4e23j,m4e23o,m4e23p,m4e23q,m4e24,m4e24d,m4h3,m4h3a,m4h4,m4h5,m4h6,m4i0m1,m4i0m2,m4i0m3,m4i0m4,m4i0m5,m4i0n1,m4i0n2,m4i0n3,m4i0n4,m4i0n5,m4i23g,m4i23i,m4j22j,m4j22k,m4k25a,m4k25b,IH5r10,IH5s2,IH5s3,IH5s4,IH5s5,o5c10,o5d2,o5d3,o5d4,o5d5,n5g1f,n5g1h,hv5_dsae,hv5_dspr,hv5_dsraw,hv5_dsss,hv5_ppvtae,hv5_ppvtpr,hv5_ppvtraw,hv5_ppvtss,p5m2a,p5m2b,p5m2c,p5m2d,p5m2e,p5m3a,p5m3b,p5m3c,p5m3d,p5m3e,p5m5a,p5m5b,p5m5c,m5a2,m5a4,m5c2b,m5c2f,m5c6a,m5c6b,m5c6c,m5c6d,m5c6e,m5c6f,m5c6g,m5c6h,m5c6i,m5c6j,m5c6o,m5c6p,m5c7,m5c8a,m5c8e,m5d2,m5d19b,m5d19f,m5d20a,m5d20b,m5d20c,m5d20d,m5d20e,m5d20f,m5d20g,m5d20h,m5d20i,m5d20j,m5d20o,m5d20p,m5d21,m5d22,m5d22d,m5e3,m5e3a,m5e4,m5e5,m5e6,m5g21a,m5g21b,m5g21c,m5g21d,m5g21e,m5g21f,m5g21g,m5g21h,m5g21i,m5g21k,m5i25a,m5i25b,p5q1c,p5q1d,p5q1f,p5q1g,p5q1h,p5q1i,p5q1j,p5q1k,p5q1m,p5q1n,p5q2a,p5q2b,p5q2c,p5q2d,p5q2e,ff_id,m1city,age3cust:age9.refusecp)
Data.PAF = FullData %>%
dplyr::select("ff_id", "k6d2b", "k6d2e", "k6d2f", "k6d2g", "k6d2h", "k6d2i", "k6d2k", "k6d2l", "k6d2m", "k6d2o", "k6d2s", "k6d2u", "k6d2v", "k6d2w", "k6d2y", "k6d2aa", "k6d2ad", "k6d2ae", "k6d2af", "k6d2ah")
Data.AnxDep = FullData %>%
dplyr::select(ff_id, k6d2ag, k6d2ai, k6d2d, k6d2j, k6d2t, k6d2ac, k6d2ak, k6d2c, k6d2n, k6d2x, p6b36, p6b40, p6b52, p6b53, p6b54, p6b68, p6b65, p6b66)
Data.SchoolConn = FullData %>%
dplyr::select("ff_id", "k5e1a", "k5e1b", "k5e1c", "k5e1d", "k6b1a", "k6b1b", "k6b1c", "k6b1d")
Extern_Variables = FullData %>%
dplyr::select(ff_id, k6d61a:k6d61m, k6d2a, k6d2p, k6d2r, k6d2z, k6d2ab, k6d2aj, k6d40, k6d48, k6f63, k6f68, k6f74, p6b35, p6b37:p6b39, p6b41:p6b45, p6b57, p6b59, p6b49:p6b51, p6b60:p6b64, p6b67)
CBCL = FullData %>%
dplyr::select(ff_id,p5q3m,p5q3ab,p5q3ac,p5q3ad,p5q3ae,p5q3af,p5q3ah,p5q3ar,p5q3av,p5q3ax,p5q3bq,p5q3ck,p5q3db,p5q3e,p5q3ao,p5q3bk,p5q3bo,p5q3bu,p5q3cu,p5q3cv,p5q3da,p5q3as,p5q3au,p5q3aw,p5q3az,p5q3bb1,p5q3bb2,p5q3bb3,p5q3bb4,p5q3bb5,p5q3bb6,p5q3bb7,p5q3b,p5q3x,p5q3aa,p5q3al,p5q3ap,p5q3bi,p5q3bm,p5q3br,p5q3bs,p5q3bz,p5q3ca,p5q3cj,p5q3cp,p5q3cr,p5q3ct,p5q3cx,p5q3cy,p5q3c,p5q3o,p5q3r,p5q3s,p5q3t,p5q3u,p5q3v,p5q3aj,p5q3bc,p5q3bn,p5q3cf,p5q3cg,p5q3ch,p5q3ci,p5q3cn,p5q3co,p5q3cq,p5q3cw)
Demo_Variables = FullData %>%
dplyr::select(ff_id, ck6ethrace, cm1bsex, cm1inpov, cm2povco, cm3povco, cm4povco, cm5povco, cp6povco)
PC_Info = FullData %>%
select(ff_id, IH3int5, IH5int5, pcg5idstat)
CBCL_15 = FullData %>%
dplyr::select(ff_id, p6b36, p6b40, p6b52, p6b53, p6b54, p6b68, p6b65, p6b66, p6b35, p6b37, p6b38, p6b39, p6b41, p6b42, p6b43, p6b44, p6b45, p6b57, p6b59, p6b49, p6b50, p6b51, p6b60, p6b61, p6b62, p6b63, p6b64, p6b67)Recode Variables
School Connectedness
We need an extra NA code for the school connectedness variables (particularly the age 15 variables) – 7 on the school connectedness scale means (NA - Homeschooled)
is.na(Data.SchoolConn) <- Data.SchoolConn == 7
Data.SchoolConn$k6b1a_r = ((Data.SchoolConn$k6b1a - 5)*-1)
Data.SchoolConn$k6b1b_r = ((Data.SchoolConn$k6b1b - 5)*-1)
Data.SchoolConn$k6b1c_r = ((Data.SchoolConn$k6b1c - 5)*-1)
Data.SchoolConn$k6b1d_r = ((Data.SchoolConn$k6b1d - 5)*-1)Teen Report Internalizing
Data.AnxDep$k6d2ag_r = ((Data.AnxDep$k6d2ag-5)*-1)
Data.AnxDep$k6d2ai_r = ((Data.AnxDep$k6d2ai-5)*-1)
Data.AnxDep$k6d2d_r = ((Data.AnxDep$k6d2d-5)*-1)
Data.AnxDep$k6d2j_r = ((Data.AnxDep$k6d2j-5)*-1)
Data.AnxDep$k6d2t_r = ((Data.AnxDep$k6d2t-5)*-1)
Data.AnxDep$k6d2ac_r = ((Data.AnxDep$k6d2ac-5)*-1)
Data.AnxDep$k6d2ak_r = ((Data.AnxDep$k6d2ak-5)*-1)
Data.AnxDep$k6d2c_r = ((Data.AnxDep$k6d2c-5)*-1)
Data.AnxDep$k6d2n_r = ((Data.AnxDep$k6d2n-5)*-1)
Data.AnxDep$k6d2x_r = ((Data.AnxDep$k6d2x-5)*-1)Teen Report Externalizing
Extern_Variables = Extern_Variables %>%
mutate(k6d2a_r = (k6d2a-5)*-1,
k6d2p_r = (k6d2p-5)*-1,
k6d2r_r = (k6d2r-5)*-1,
k6d2z_r = (k6d2z-5)*-1,
k6d2ab_r = (k6d2ab-5)*-1,
k6d2aj_r = (k6d2aj-5)*-1,
k6d40_r = (k6d40-3)*-1,
k6d48_r = (k6d48-3)*-1,
k6f63_r = (k6f63-3)*-1,
k6f68_r = (k6f68-3)*-1,
k6f74_r = (k6f74-2)*-1)Positive Adolescent Function
Data.PAF$k6d2b_r = ((Data.PAF$k6d2b-5)*-1)
Data.PAF$k6d2b_r = 5-Data.PAF$k6d2b
Data.PAF$k6d2e_r = ((Data.PAF$k6d2e-5)*-1)
Data.PAF$k6d2f_r = ((Data.PAF$k6d2f-5)*-1)
Data.PAF$k6d2g_r = ((Data.PAF$k6d2g-5)*-1)
Data.PAF$k6d2h_r = ((Data.PAF$k6d2h-5)*-1)
Data.PAF$k6d2i_r = ((Data.PAF$k6d2i-5)*-1)
Data.PAF$k6d2k_r = ((Data.PAF$k6d2k-5)*-1)
Data.PAF$k6d2l_r = ((Data.PAF$k6d2l-5)*-1)
Data.PAF$k6d2m_r = ((Data.PAF$k6d2m-5)*-1)
Data.PAF$k6d2o_r = ((Data.PAF$k6d2o-5)*-1)
Data.PAF$k6d2s_r = ((Data.PAF$k6d2s-5)*-1)
Data.PAF$k6d2u_r = ((Data.PAF$k6d2u-5)*-1)
Data.PAF$k6d2v_r = ((Data.PAF$k6d2v-5)*-1)
Data.PAF$k6d2w_r = ((Data.PAF$k6d2w-5)*-1)
Data.PAF$k6d2y_r = ((Data.PAF$k6d2y-5)*-1)
Data.PAF$k6d2aa_r = ((Data.PAF$k6d2aa-5)*-1)
Data.PAF$k6d2ad_r = ((Data.PAF$k6d2ad-5)*-1)
Data.PAF$k6d2ae_r = ((Data.PAF$k6d2ae-5)*-1)
Data.PAF$k6d2af_r = ((Data.PAF$k6d2af-5)*-1)
Data.PAF$k6d2ah_r = ((Data.PAF$k6d2ah-5)*-1)Demographic Variable Prep
# Race/Ethnicity at Age 15- ck6ethrace
# 5 Multi-racial, non-hispanic
# 4 Other only, non-hispanic
# 3 Hispanic/Latino
# 2 Black/Af. American only, non-hispanic
# 1 White only, non-hispanic
Demo_Variables = Demo_Variables %>%
mutate(povco_avg = (cm1inpov + cm2povco + cm3povco + cm4povco + cm5povco + cp6povco)/6) %>%
mutate(Race_AA = if_else(ck6ethrace == 2, 1, 0),
Race_C = if_else(ck6ethrace == 1, 1, 0),
Race_L = if_else(ck6ethrace == 3, 1, 0)) %>%
mutate(cm1bsex = cm1bsex-1) %>%
dplyr::select(ff_id, povco_avg, Race_AA, Race_C, Race_L, cm1bsex)Merge to One Dataframe
All_Variables = Data.ThreatDepVars %>%
left_join(Data.AnxDep, by="ff_id")
All_Variables = All_Variables %>%
left_join(Extern_Variables, by="ff_id")
All_Variables = All_Variables %>%
left_join(Data.PAF, by = "ff_id")
All_Variables = All_Variables %>%
left_join(Data.SchoolConn, by = "ff_id")
All_Variables = All_Variables %>%
left_join(Demo_Variables, by = "ff_id")Write Variables to CSV
write.csv(All_Variables,'/Volumes/csmonk/csmonk-lab/Data/Fragile_Families_Longitudinal/Leigh_FFCWS/VE_SD_IntExtPAF_SchoolConn_090420.csv')Write Variables for MPlus
All_Variables[is.na(All_Variables)] = 99prepareMplusData(All_Variables,'../../Manuscript/Data_Analysis/All_Variables_090420.dat')CBCL - Age 9
V.M. suggested that we control for psychopathology at Age 9 since we have the data. We don’t have a teen report of psychopathology at this age, but we do have parent report in the CBCL. To not make the SEM too complex, I am just feeding in sum scores representing internalizing and externalizing.
CBCL_AnxDep=p5q3m,p5q3ab,p5q3ac,p5q3ad,p5q3ae,p5q3af,p5q3ah,p5q3ar,p5q3av,p5q3ax,p5q3bq,p5q3ck,p5q3db
CBCL_With=p5q3e,p5q3ao,p5q3bk,p5q3bo,p5q3bu,p5q3cu,p5q3cv,p5q3da
CBCL_Som=p5q3as,p5q3au,p5q3aw,p5q3az,p5q3bb1,p5q3bb2,p5q3bb3,p5q3bb4,p5q3bb5,p5q3bb6,p5q3bb7
CBCL_Rule=p5q3b,p5q3x,p5q3aa,p5q3al,p5q3ap,p5q3bi,p5q3bm,p5q3br,p5q3bs,p5q3bz,p5q3ca,p5q3cj,p5q3cp,p5q3cr,p5q3ct,p5q3cx,p5q3cy
CBCL_Agg=p5q3c,p5q3o,p5q3r,p5q3s,p5q3t,p5q3u,p5q3v,p5q3aj,p5q3bc,p5q3bn,p5q3cf,p5q3cg,p5q3ch,p5q3ci,p5q3cn,p5q3co,p5q3cq,p5q3cw
CBCL = CBCL %>%
mutate(InternCBCL9 = p5q3m + p5q3ab + p5q3ac + p5q3ad + p5q3ae + p5q3af + p5q3ah + p5q3ar + p5q3av + p5q3ax + p5q3bq + p5q3ck + p5q3db + p5q3e + p5q3ao + p5q3bk + p5q3bo + p5q3bu + p5q3cu + p5q3cv + p5q3da + p5q3as + p5q3au + p5q3aw + p5q3az + p5q3bb1 + p5q3bb2 + p5q3bb3 + p5q3bb4 + p5q3bb5 + p5q3bb6 + p5q3bb7,
ExternCBCL9 = p5q3b + p5q3x + p5q3aa + p5q3al + p5q3ap + p5q3bi + p5q3bm + p5q3br + p5q3bs + p5q3bz + p5q3ca + p5q3cj + p5q3cp + p5q3cr + p5q3ct + p5q3cx + p5q3cy + p5q3c + p5q3o + p5q3r + p5q3s + p5q3t + p5q3u + p5q3v + p5q3aj + p5q3bc + p5q3bn + p5q3cf + p5q3cg + p5q3ch + p5q3ci + p5q3cn + p5q3co + p5q3cq + p5q3cw)CBCL = CBCL %>%
select(ff_id, InternCBCL9, ExternCBCL9)
All_Variables = All_Variables %>%
left_join(CBCL, by='ff_id')CBCL - Age 15
CBCL_15 = CBCL_15 %>%
mutate(Intern_CBCL15 = p6b36 + p6b40 + p6b52 + p6b53 + p6b54 + p6b68 + p6b65 + p6b66,
Extern_CBCL15 = p6b35 + p6b37 + p6b38 + p6b39 + p6b41 + p6b42 + p6b43 + p6b44 + p6b45 + p6b57 + p6b59 + p6b49 + p6b50 + p6b51 + p6b60 + p6b61 + p6b62 + p6b63 + p6b64 + p6b67) %>%
dplyr::select(ff_id, Intern_CBCL15, Extern_CBCL15)
All_Variables = All_Variables %>%
left_join(CBCL_15, by='ff_id')I wrote this out as a separate Mplus .dat
All_Variables[is.na(All_Variables)] = 99prepareMplusData(All_Variables,'../../Manuscript/Data_Analysis/All_Variables_101220_wCBCL.dat')
prepareMplusData(All_Variables,'../../Manuscript/Data_Analysis/All_Variables_101320_wCBCL_forSC9Check.dat')Actual Analysis Starts Here - Read in Data
Here I am reading in data written out from the creation/recoding of variables above as well as the factor scores extracted from Mplus for the simple slopes analysis. Convention for how factor scores should be extracted when probing latent variable interactions has not been clearly established. Therefore, we used two different approaches and tested whether the method used influenced the interpretation of the interaction. The first method extracted the factor scores from a measurement model containing all of the latent variables in the model (e.g. school connectedness at age 9, school connectedness at age 15, and positive adolescent function). The second approach used the factor scores extracted from individual measurement models. Due to the correlations between latent variables that are accounted for in measurement models (Kline, 2015), it is likely that the extracted factors are different depending on how they are extracted. This is why there are two dataframes, All_Variables and All_Variables2.
# OLD VERSION _ All_Variables = read_csv('/Volumes/csmonk/csmonk-lab/Data/Fragile_Families_Longitudinal/Leigh_FFCWS/VE_SD_IntExtPAF_SchoolConn_090420.csv')
All_Variables = read_csv('/Volumes/csmonk/csmonk-lab/Data/Fragile_Families_Longitudinal/Leigh_FFCWS/VE_SD_IntExtPAF_SchoolConn_CBCL915_101220.csv')## Parsed with column specification:
## cols(
## .default = col_double()
## )## See spec(...) for full column specifications.fs_scores = read_csv('FactorScores_SC159PAF_090320.csv')## Parsed with column specification:
## cols(
## ff_id = col_character(),
## SC15 = col_double(),
## SC15_SE = col_double(),
## SC9 = col_double(),
## SC9_SE = col_double(),
## PAF = col_double(),
## PAF_SE = col_double()
## )fs_scores$ff_id = as.integer(fs_scores$ff_id)## Warning: NAs introduced by coercionSC9 = read_csv('SC9_080520.csv')## Parsed with column specification:
## cols(
## ff_id = col_character(),
## SC9 = col_double(),
## SC9_SE = col_double()
## )SC9$ff_id = as.integer(SC9$ff_id)## Warning: NAs introduced by coercionSC15 = read_csv('SC15_080520.csv')## Parsed with column specification:
## cols(
## ff_id = col_double(),
## SC15 = col_double(),
## SC15_SE = col_double()
## )SC15$ff_id = as.integer(SC15$ff_id)
Int9 = read_csv('Int9_102720.csv')## Parsed with column specification:
## cols(
## ff_id = col_character(),
## Int9 = col_double(),
## Int9_SE = col_double()
## )Int9$ff_id = as.integer(Int9$ff_id)## Warning: NAs introduced by coercionExt9 = read_csv('Ext9_102720.csv')## Parsed with column specification:
## cols(
## ff_id = col_character(),
## Ext9 = col_double(),
## Ext9_SE = col_double()
## )Ext9$ff_id = as.integer(Ext9$ff_id)## Warning: NAs introduced by coercionInt15 = read_csv('Int_080520.csv')## Parsed with column specification:
## cols(
## ff_id = col_double(),
## Int = col_double(),
## Int_SE = col_double()
## )Int15$ff_id = as.integer(Int15$ff_id)
Ext15 = read_csv('Ext_080520.csv')## Parsed with column specification:
## cols(
## ff_id = col_double(),
## Ext = col_double(),
## Ext_SE = col_double()
## )Ext15$ff_id = as.integer(Ext15$ff_id)
Ext15_red = read_csv('Ext15_reduced_102720.csv')## Parsed with column specification:
## cols(
## ff_id = col_double(),
## Ext15 = col_double(),
## Ext15_SE = col_double()
## )Ext15_red$ff_id = as.integer(Ext15_red$ff_id)
PAF = read_csv('PAF_080520.csv')## Parsed with column specification:
## cols(
## ff_id = col_double(),
## PAF = col_double(),
## PAF_SE = col_double()
## )PAF$ff_id = as.integer(PAF$ff_id)
# Merge things
fs_scores_indiv = All_Variables %>%
select(ff_id)
fs_scores_indiv = fs_scores_indiv %>%
left_join(SC9, by="ff_id")
fs_scores_indiv = fs_scores_indiv %>%
left_join(SC15, by="ff_id")
fs_scores_indiv = fs_scores_indiv %>%
left_join(Int9, by="ff_id")
fs_scores_indiv = fs_scores_indiv %>%
left_join(Ext9, by="ff_id")
fs_scores_indiv = fs_scores_indiv %>%
left_join(Int15, by="ff_id")
fs_scores_indiv = fs_scores_indiv %>%
left_join(Ext15, by="ff_id")
fs_scores_indiv = fs_scores_indiv %>%
left_join(Ext15_red, by="ff_id")
fs_scores_indiv = fs_scores_indiv %>%
left_join(PAF, by="ff_id")
# It is important that all_variables2 is created first so that there are not doubles of the factor scores.
All_Variables2 = All_Variables %>%
left_join(fs_scores_indiv, by="ff_id")
# Join everything together in the All_Variables dataframe.
All_Variables = All_Variables %>%
left_join(fs_scores, by="ff_id")
# This is adding the internalizing and externalizing factor scores separately into the All_Variables dataframe because they're not in the significant moderation model that we probe with the simple slopes but I want to show the mean/sd of the factor scores.
All_Variables = All_Variables %>%
left_join(Int9, by="ff_id") %>%
left_join(Ext9, by="ff_id") %>%
left_join(Int15, by="ff_id") %>%
left_join(Ext15_red, by="ff_id") %>%
left_join(Ext15, by="ff_id")This reverses any 99s that are marked in the MPlus data because MPlus doesn’t take NAs.
is.na(All_Variables) <- All_Variables == 99
is.na(All_Variables2) <- All_Variables2 == 99 Look at distributions of our variables
library(ggpubr)
PAF = All_Variables %>%
ggplot(mapping = aes(x=PAF)) + xlab("Positive Adolescent Function") + geom_density() + theme_classic()
Dep = All_Variables %>%
ggplot(mapping = aes(x=DepComp)) + xlab("Social Deprivation") + geom_density() + theme_classic()
Threat = All_Variables %>%
ggplot(mapping = aes(x=ThreatComp)) + xlab("Violence Exposure") + geom_density() + theme_classic()
SC9 = All_Variables %>%
ggplot(mapping = aes(x=SC9)) + xlab("School Connectedness Age 9") + geom_density() + theme_classic()
SC15 = All_Variables %>%
ggplot(mapping = aes(x=SC15)) + xlab("School Connectedness Age 15") + geom_density() + theme_classic()
Int9 = All_Variables %>%
ggplot(mapping = aes(x=Int9)) + xlab("Internalizing Psychopathology Age 9") + geom_density() + theme_classic()
Ext9 = All_Variables %>%
ggplot(mapping = aes(x=Ext9)) + xlab("Externalizing Psychopathology Age 9") + geom_density() + theme_classic()
Int15 = All_Variables %>%
ggplot(mapping = aes(x=Int)) + xlab("Internalizing Psychopathology Age 15") + geom_density() + theme_classic()
Ext15 = All_Variables %>%
ggplot(mapping = aes(x=Ext15)) + xlab("Externalizing Psychopathology Age 15") + geom_density() + theme_classic()
figure <- ggarrange(PAF, Threat, Dep, SC9,SC15, Int9, Ext9, Int15,Ext15,
ncol = 2, nrow = 6)## Warning: Removed 1168 rows containing non-finite values (stat_density).## Warning: Removed 382 rows containing non-finite values (stat_density).## Warning: Removed 337 rows containing non-finite values (stat_density).## Warning: Removed 1168 rows containing non-finite values (stat_density).
## Warning: Removed 1168 rows containing non-finite values (stat_density).## Warning: Removed 1562 rows containing non-finite values (stat_density).
## Warning: Removed 1562 rows containing non-finite values (stat_density).## Warning: Removed 1461 rows containing non-finite values (stat_density).## Warning: Removed 1460 rows containing non-finite values (stat_density).figure
Skewness and Kurtosis
library(psych)##
## Attaching package: 'psych'## The following objects are masked from 'package:ggplot2':
##
## %+%, alphapsych::describe(All_Variables$ThreatComp) # -- this is potentially problematic, it seems like there are some serious kurtosis issues ## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 4516 0.01 0.63 -0.13 -0.07 0.43 -1.2 13.26 14.46 4.55 59.84 0.01psych::describe(All_Variables$DepComp)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 4561 0.01 0.57 -0.08 -0.04 0.49 -1.58 4.02 5.6 1.19 3.17 0.01psych::describe(All_Variables$PAF)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 3730 -0.03 0.85 -0.02 -0.02 0.83 -3.08 1.7 4.78 -0.15 -0.12 0.01psych::describe(All_Variables$Int9)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 3336 0.09 0.82 0.05 0.03 0.85 -1.01 4.91 5.92 0.68 1.15 0.01psych::describe(All_Variables$Ext9)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 3336 0.06 0.88 0.03 0.01 0.93 -1.2 4.28 5.48 0.44 0.02 0.02psych::describe(All_Variables$Int)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 3437 0.03 0.88 0.02 0.01 0.95 -1.38 3.19 4.57 0.24 -0.46 0.01psych::describe(All_Variables$Ext15)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 3438 0.04 0.9 0.07 0.06 0.89 -2.05 2.76 4.8 -0.15 -0.2 0.02psych::describe(All_Variables$SC9)## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 3730 -0.05 0.74 0 0.01 0.74 -2.43 1.02 3.45 -0.64 0.14
## se
## X1 0.01psych::describe(All_Variables$SC15)## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 3730 -0.04 0.78 -0.01 0 0.78 -2.88 1.34 4.22 -0.46 0
## se
## X1 0.01Table Creation
The code below creates: - Demographic Table - Correlation Table w/descriptives for our main variables.
# Import a categorical race variable because I had accidentally cut it out earlier.
Race_Cat = read_csv('Race_Categorical.csv')## Parsed with column specification:
## cols(
## ff_id = col_double(),
## Race_Cat = col_character()
## )Race_Cat$Race_Cat = as.factor(Race_Cat$Race_Cat)
MatBirthVars = read_csv('Maternal_Birth_Variables.csv')## Parsed with column specification:
## cols(
## ff_id = col_double(),
## m1b2 = col_double(),
## cm1edu = col_double()
## )MatBirthVars$m1b2 = as.factor(MatBirthVars$m1b2)
MatBirthVars$cm1edu = as.factor(MatBirthVars$cm1edu)
All_Variables = All_Variables %>%
left_join(Race_Cat, by='ff_id') %>%
left_join(MatBirthVars, by='ff_id')
All_Variables2 = All_Variables2 %>%
left_join(Race_Cat, by='ff_id') %>%
left_join(MatBirthVars, by='ff_id')
# This sex variable is created for the demographic table
All_Variables = All_Variables %>%
mutate(sex = case_when(
cm1bsex == 0 ~ 'Male',
cm1bsex == 1 ~ 'Female',
TRUE ~ 'NA'
))
All_Variables = All_Variables %>%
mutate(MatMar = case_when(
m1b2 == 1 ~ 'Married',
m1b2 == 2 ~ 'Not Married',
TRUE ~ 'NA'
),
MatEdu = case_when(
cm1edu == 1 ~ 'Less than high school',
cm1edu == 2 ~ 'High school or equivalent',
cm1edu == 3 ~ 'Some college or technical school',
cm1edu == 4 ~ 'College or graduate school',
TRUE ~ 'NA'
))
All_Variables2 = All_Variables2 %>%
mutate(sex = case_when(
cm1bsex == 0 ~ 'Male',
cm1bsex == 1 ~ 'Female',
TRUE ~ 'NA'
),
MatEdu = case_when(
cm1edu == 1 ~ 'Less than high school',
cm1edu == 2 ~ 'High school or equivalent',
cm1edu == 3 ~ 'Some college or technical school',
cm1edu == 4 ~ 'College or graduate school',
TRUE ~ 'NA'
))library(table1)##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<-library(MatchIt)
table1::label(All_Variables$sex) = "Sex"
table1::label(All_Variables$Race_Cat) = "Race-Ethnicity"
table1::label(All_Variables$povco_avg) = "Average Poverty Ratio"
table1::label(All_Variables$MatMar) = "Maternal Marital Status at Child's Birth"
table1::label(All_Variables$MatEdu) = "Maternal Education at Child's Birth"
table1::label(All_Variables$ThreatComp) = "Violence Exposure"
table1::label(All_Variables$DepComp) = "Social Deprivation"
table1::label(All_Variables$SC9) = "School Connectedness Age 9"
table1::label(All_Variables$SC15) = "School Connectedness Age 15"
table1::label(All_Variables$PAF) = "Positive Adolescent Function"
table1::label(All_Variables$Int) = "Internalizing Symptoms Age 15"
table1::label(All_Variables$Ext15) = "Externalizing Symptoms Age 15"
table1::label(All_Variables$Int9) = "Internalizing Symptoms Age 9"
table1::label(All_Variables$Ext9) = "Externalizing Symptoms Age 9"
# This next piece is trying to use the people that would be in the MPlus model only
fs_scores_id = fs_scores %>%
select("ff_id")
All_Variables_Red = All_Variables %>%
right_join(fs_scores_id, by='ff_id')
table1::table1(~sex + Race_Cat + povco_avg + MatMar + MatEdu, data = All_Variables_Red)| Overall (N=3731) |
|
|---|---|
| Sex | |
| Female | 1787 (47.9%) |
| Male | 1944 (52.1%) |
| Race-Ethnicity | |
| African American | 1601 (42.9%) |
| Caucasian | 590 (15.8%) |
| Latinx | 813 (21.8%) |
| Other | 261 (7.0%) |
| Missing | 466 (12.5%) |
| Average Poverty Ratio | |
| Mean (SD) | 2.12 (2.06) |
| Median [Min, Max] | 1.48 [0.121, 21.2] |
| Missing | 995 (26.7%) |
| Maternal Marital Status at Child's Birth | |
| Married | 881 (23.6%) |
| NA | 21 (0.6%) |
| Not Married | 2829 (75.8%) |
| Maternal Education at Child's Birth | |
| College or graduate school | 399 (10.7%) |
| High school or equivalent | 1178 (31.6%) |
| Less than high school | 1218 (32.6%) |
| NA | 6 (0.2%) |
| Some college or technical school | 930 (24.9%) |
table1::table1(~ThreatComp + DepComp + SC9 + SC15 + Int9 + Ext9 + Int + Ext15 + PAF, data = All_Variables)| Overall (N=4898) |
|
|---|---|
| Violence Exposure | |
| Mean (SD) | 0.00513 (0.632) |
| Median [Min, Max] | -0.132 [-1.20, 13.3] |
| Missing | 382 (7.8%) |
| Social Deprivation | |
| Mean (SD) | 0.00522 (0.568) |
| Median [Min, Max] | -0.0833 [-1.58, 4.02] |
| Missing | 337 (6.9%) |
| School Connectedness Age 9 | |
| Mean (SD) | -0.0462 (0.735) |
| Median [Min, Max] | 0.00400 [-2.43, 1.02] |
| Missing | 1168 (23.8%) |
| School Connectedness Age 15 | |
| Mean (SD) | -0.0409 (0.777) |
| Median [Min, Max] | -0.00600 [-2.88, 1.34] |
| Missing | 1168 (23.8%) |
| Internalizing Symptoms Age 9 | |
| Mean (SD) | 0.0898 (0.820) |
| Median [Min, Max] | 0.0490 [-1.01, 4.91] |
| Missing | 1562 (31.9%) |
| Externalizing Symptoms Age 9 | |
| Mean (SD) | 0.0589 (0.875) |
| Median [Min, Max] | 0.0330 [-1.20, 4.28] |
| Missing | 1562 (31.9%) |
| Internalizing Symptoms Age 15 | |
| Mean (SD) | 0.0288 (0.877) |
| Median [Min, Max] | 0.0200 [-1.38, 3.19] |
| Missing | 1461 (29.8%) |
| Externalizing Symptoms Age 15 | |
| Mean (SD) | 0.0414 (0.898) |
| Median [Min, Max] | 0.0705 [-2.05, 2.76] |
| Missing | 1460 (29.8%) |
| Positive Adolescent Function | |
| Mean (SD) | -0.0260 (0.855) |
| Median [Min, Max] | -0.0160 [-3.08, 1.70] |
| Missing | 1168 (23.8%) |
CorTable = All_Variables %>%
dplyr::select(ThreatComp, DepComp, SC9, SC15, Int9, Ext9, Int, Ext15, PAF)
library(apaTables)
apa.cor.table(CorTable, filename = '/Users/leighgoetschius/Box Sync/Shift_Persist_FF/Manuscript/Figures/cor_table_102720.doc', table.number = NA,
show.conf.interval = TRUE, landscape = TRUE)Creating some interaction terms
All_Variables = All_Variables %>%
mutate(DepSC9 = SC9*DepComp,
DepSC15 = SC15*DepComp,
ThrSC9 = SC9*ThreatComp,
ThrSC15 = SC15*ThreatComp)
All_Variables2 = All_Variables2 %>%
mutate(DepSC9 = SC9*DepComp,
DepSC15 = SC15*DepComp,
ThrSC9 = SC9*ThreatComp,
ThrSC15 = SC15*ThreatComp)Looking at the SEM model in regression space
Using factor scores extracted from the group CFA
summary(lm(PAF ~ DepComp + ThreatComp, data = All_Variables))##
## Call:
## lm(formula = PAF ~ DepComp + ThreatComp, data = All_Variables)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.05025 -0.53320 0.00388 0.57041 1.92635
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.02442 0.01394 -1.752 0.0799 .
## DepComp -0.17690 0.02866 -6.172 7.48e-10 ***
## ThreatComp -0.02166 0.02826 -0.766 0.4436
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8493 on 3708 degrees of freedom
## (1187 observations deleted due to missingness)
## Multiple R-squared: 0.01371, Adjusted R-squared: 0.01318
## F-statistic: 25.78 on 2 and 3708 DF, p-value: 7.628e-12summary(lm(PAF ~ DepComp + ThreatComp + SC9 + SC15, data = All_Variables))##
## Call:
## lm(formula = PAF ~ DepComp + ThreatComp + SC9 + SC15, data = All_Variables)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7386 -0.4077 -0.0037 0.4152 2.3248
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.004281 0.010519 0.407 0.68407
## DepComp -0.070575 0.021666 -3.257 0.00113 **
## ThreatComp 0.067210 0.021336 3.150 0.00165 **
## SC9 0.024069 0.014924 1.613 0.10689
## SC15 0.721074 0.014277 50.506 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6392 on 3706 degrees of freedom
## (1187 observations deleted due to missingness)
## Multiple R-squared: 0.4417, Adjusted R-squared: 0.4411
## F-statistic: 732.9 on 4 and 3706 DF, p-value: < 2.2e-16summary(lm(PAF ~ DepComp + SC9 + DepSC9 + ThreatComp + SC15, data = All_Variables))##
## Call:
## lm(formula = PAF ~ DepComp + SC9 + DepSC9 + ThreatComp + SC15,
## data = All_Variables)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.74000 -0.40724 -0.00091 0.41538 2.34927
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.002402 0.010534 0.228 0.819646
## DepComp -0.074880 0.021708 -3.449 0.000568 ***
## SC9 0.024531 0.014913 1.645 0.100071
## DepSC9 -0.074438 0.027981 -2.660 0.007840 **
## ThreatComp 0.067758 0.021319 3.178 0.001494 **
## SC15 0.721054 0.014265 50.546 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6386 on 3705 degrees of freedom
## (1187 observations deleted due to missingness)
## Multiple R-squared: 0.4427, Adjusted R-squared: 0.442
## F-statistic: 588.7 on 5 and 3705 DF, p-value: < 2.2e-16summary(lm(PAF ~ DepComp + SC9 + SC15 + DepSC9 + DepSC15, data = All_Variables))##
## Call:
## lm(formula = PAF ~ DepComp + SC9 + SC15 + DepSC9 + DepSC15, data = All_Variables)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.73527 -0.40654 -0.00112 0.41710 2.28553
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0007214 0.0106254 0.068 0.94587
## DepComp -0.0569283 0.0199981 -2.847 0.00444 **
## SC9 0.0238164 0.0149430 1.594 0.11106
## SC15 0.7166728 0.0142452 50.310 < 2e-16 ***
## DepSC9 -0.0529355 0.0295423 -1.792 0.07324 .
## DepSC15 -0.0594081 0.0265592 -2.237 0.02536 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6399 on 3708 degrees of freedom
## (1184 observations deleted due to missingness)
## Multiple R-squared: 0.441, Adjusted R-squared: 0.4402
## F-statistic: 585 on 5 and 3708 DF, p-value: < 2.2e-16Using factor scores extracted from the individual CFAs
summary(lm(PAF ~ DepComp + ThreatComp, data = All_Variables2))##
## Call:
## lm(formula = PAF ~ DepComp + ThreatComp, data = All_Variables2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.03521 -0.59759 -0.01079 0.64018 1.81515
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.02167 0.01512 -1.434 0.152
## DepComp -0.18102 0.03087 -5.864 4.95e-09 ***
## ThreatComp -0.00880 0.03064 -0.287 0.774
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8836 on 3415 degrees of freedom
## (1480 observations deleted due to missingness)
## Multiple R-squared: 0.01242, Adjusted R-squared: 0.01184
## F-statistic: 21.48 on 2 and 3415 DF, p-value: 5.39e-10summary(lm(PAF ~ DepComp + ThreatComp + SC9 + SC15, data = All_Variables2))##
## Call:
## lm(formula = PAF ~ DepComp + ThreatComp + SC9 + SC15, data = All_Variables2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.57499 -0.55960 -0.03053 0.54348 2.58996
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.004298 0.014682 0.293 0.769726
## DepComp -0.110800 0.030553 -3.627 0.000292 ***
## ThreatComp 0.058713 0.031050 1.891 0.058730 .
## SC9 0.061349 0.019283 3.182 0.001480 **
## SC15 0.456896 0.019134 23.879 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8005 on 2989 degrees of freedom
## (1904 observations deleted due to missingness)
## Multiple R-squared: 0.1812, Adjusted R-squared: 0.1801
## F-statistic: 165.4 on 4 and 2989 DF, p-value: < 2.2e-16summary(lm(PAF ~ DepComp + SC9 + DepSC9 + ThreatComp + SC15, data = All_Variables2))##
## Call:
## lm(formula = PAF ~ DepComp + SC9 + DepSC9 + ThreatComp + SC15,
## data = All_Variables2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.52217 -0.56513 -0.02159 0.55240 2.62476
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.002312 0.014693 0.157 0.87499
## DepComp -0.117474 0.030655 -3.832 0.00013 ***
## SC9 0.062400 0.019273 3.238 0.00122 **
## DepSC9 -0.088518 0.036923 -2.397 0.01657 *
## ThreatComp 0.058182 0.031026 1.875 0.06085 .
## SC15 0.457154 0.019119 23.911 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7999 on 2988 degrees of freedom
## (1904 observations deleted due to missingness)
## Multiple R-squared: 0.1828, Adjusted R-squared: 0.1814
## F-statistic: 133.7 on 5 and 2988 DF, p-value: < 2.2e-16summary(lm(PAF ~ DepComp + ThreatComp + SC9 + SC15 + DepSC9 + DepSC15, data = All_Variables2))##
## Call:
## lm(formula = PAF ~ DepComp + ThreatComp + SC9 + SC15 + DepSC9 +
## DepSC15, data = All_Variables2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.56640 -0.56484 -0.02891 0.55590 2.56395
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.001603 0.014806 -0.108 0.91377
## DepComp -0.127004 0.030979 -4.100 4.25e-05 ***
## ThreatComp 0.058162 0.031009 1.876 0.06080 .
## SC9 0.062309 0.019262 3.235 0.00123 **
## SC15 0.458755 0.019124 23.989 < 2e-16 ***
## DepSC9 -0.074683 0.037497 -1.992 0.04649 *
## DepSC15 -0.073970 0.035568 -2.080 0.03764 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7995 on 2987 degrees of freedom
## (1904 observations deleted due to missingness)
## Multiple R-squared: 0.184, Adjusted R-squared: 0.1824
## F-statistic: 112.3 on 6 and 2987 DF, p-value: < 2.2e-16Plot of the SC15 Effect
PAF_form = lm(All_Variables2$PAF~All_Variables2$SC15)
Int_form = lm(All_Variables2$Int~All_Variables2$SC15)
Ext_form = lm(All_Variables2$Ext15~All_Variables2$SC15)
PAF_SC15 = All_Variables2 %>%
ggplot(aes(x=SC15, y=PAF)) +
geom_point(position="jitter") +
geom_smooth(method="lm", colour="skyblue3") +
xlab("School Connectedness Age 15") +
ylab("Positive Adolescent Function Age 15") +
stat_regline_equation(aes(label = paste(..adj.rr.label..)),formula = PAF_form) +
theme_classic() +
theme(text = element_text(family='serif'))
Int_SC15 = All_Variables2 %>%
ggplot(aes(x=SC15, y=Int)) +
geom_point(position="jitter") +
geom_smooth(method="lm", colour="red") +
xlab("School Connectedness Age 15") +
ylab("Internalizing Symptoms Age 15") +
stat_regline_equation(aes(label = paste(..adj.rr.label..)),formula = Int_form) +
theme_classic() +
theme(text = element_text(family='serif'))
Ext_SC15 =All_Variables2 %>%
ggplot(aes(x=SC15, y=Ext15)) +
geom_point(position="jitter") +
geom_smooth(method="lm", colour = "red")+
xlab("School Connectedness Age 15") +
ylab("Externalizing Symptoms Age 15") +
stat_regline_equation(aes(label = paste(..adj.rr.label..)),formula = Ext_form) +
theme_classic() +
theme(text = element_text(family='serif'))
SC15_figure <- ggarrange(Int_SC15, Ext_SC15, PAF_SC15,
ncol = 2, nrow = 2)## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 1518 rows containing non-finite values (stat_smooth).## Warning: Removed 1518 rows containing non-finite values (stat_regline_equation).## Warning: Removed 1518 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 1517 rows containing non-finite values (stat_smooth).## Warning: Removed 1517 rows containing non-finite values (stat_regline_equation).## Warning: Removed 1517 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 1518 rows containing non-finite values (stat_smooth).## Warning: Removed 1518 rows containing non-finite values (stat_regline_equation).## Warning: Removed 1518 rows containing missing values (geom_point).SC15_figure
Plot SC9 Contemporaneous Effect
Int9_form = lm(All_Variables2$Int9~All_Variables2$SC9)
Ext9_form = lm(All_Variables2$Ext9~All_Variables2$SC9)
Int9_SC9 = All_Variables2 %>%
ggplot(aes(x=SC9, y=Int9)) +
geom_point(position="jitter") +
geom_smooth(method="lm", colour="red") +
xlab("School Connectedness Age 9") +
ylab("Internalizing Symptoms Age 9") +
stat_regline_equation(aes(label = paste(..adj.rr.label..)),formula = Int9_form) +
theme_classic() +
theme(text = element_text(family='serif'))
Ext9_SC9 =All_Variables2 %>%
ggplot(aes(x=SC9, y=Ext9)) +
geom_point(position="jitter") +
geom_smooth(method="lm", colour = "red")+
xlab("School Connectedness Age 9") +
ylab("Externalizing Symptoms Age 9") +
stat_regline_equation(aes(label = paste(..adj.rr.label..)),formula = Ext9_form) +
theme_classic() +
theme(text = element_text(family='serif'))
SC9_figure <- ggarrange(Int9_SC9, Ext9_SC9,
ncol = 2, nrow = 1)## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 1626 rows containing non-finite values (stat_smooth).## Warning: Removed 1626 rows containing non-finite values (stat_regline_equation).## Warning: Removed 1626 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 1626 rows containing non-finite values (stat_smooth).## Warning: Removed 1626 rows containing non-finite values (stat_regline_equation).## Warning: Removed 1626 rows containing missing values (geom_point).SC9_figure
Probing the SDxSC9 Interaction
# Center/Standardize the Variables
All_Variables = All_Variables %>%
mutate(PAF_scaled = scale(PAF),
DepComp_scaled = scale(DepComp),
SC9_scaled = scale(SC9))
# Set your model to a variable
interact_mod = lm(PAF_scaled ~ DepComp_scaled * SC9_scaled, data = All_Variables)
# Interpret the interaction
interact_plot(interact_mod, pred = DepComp_scaled, modx = SC9_scaled, interval = FALSE, rug = TRUE, rug.sides = 'bl', modx.labels = c("Low","Mean","High"), x.label = "Social Deprivation", y.label = "Positive Adolescent Function", legend.main = "School Connectedness", data = All_Variables) + geom_rect(aes(xmin=-2.79, xmax = 2.83, ymin=-4, ymax = 2), fill = 'light gray', alpha = 0.005, color = 'black') + theme_classic() + theme(text=element_text(size = 16, family="serif"))
johnson_neyman(interact_mod, pred = DepComp_scaled, modx = SC9_scaled)## JOHNSON-NEYMAN INTERVAL
##
## When SC9_scaled is OUTSIDE the interval [-41.78, -1.41], the slope of
## DepComp_scaled is p < .05.
##
## Note: The range of observed values of SC9_scaled is [-3.25, 1.45]
johnson_neyman(interact_mod, pred = SC9_scaled, modx = DepComp_scaled)## JOHNSON-NEYMAN INTERVAL
##
## When DepComp_scaled is OUTSIDE the interval [2.83, 78.00], the slope of
## SC9_scaled is p < .05.
##
## Note: The range of observed values of DepComp_scaled is [-2.79, 7.07]
sim_slopes(interact_mod, pred = DepComp_scaled, modx = SC9_scaled, data = All_Variables, digits = 3)## JOHNSON-NEYMAN INTERVAL
##
## When SC9_scaled is OUTSIDE the interval [-41.776, -1.415], the slope of
## DepComp_scaled is p < .05.
##
## Note: The range of observed values of SC9_scaled is [-3.246, 1.445]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of DepComp_scaled when SC9_scaled = -1.002 (- 1 SD):
##
## Est. S.E. t val. p
## -------- ------- -------- -------
## -0.075 0.024 -3.053 0.002
##
## Slope of DepComp_scaled when SC9_scaled = 0.000 (Mean):
##
## Est. S.E. t val. p
## -------- ------- -------- -------
## -0.112 0.017 -6.628 0.000
##
## Slope of DepComp_scaled when SC9_scaled = 1.002 (+ 1 SD):
##
## Est. S.E. t val. p
## -------- ------- -------- -------
## -0.150 0.025 -6.065 0.000# Getting coefficients, variances, and covariances for the Preacher interaction tool.
vcov(interact_mod)## (Intercept) DepComp_scaled SC9_scaled
## (Intercept) 2.552404e-04 1.744482e-07 -3.823368e-08
## DepComp_scaled 1.744482e-07 2.865047e-04 1.754375e-05
## SC9_scaled -3.823368e-08 1.754375e-05 2.541336e-04
## DepComp_scaled:SC9_scaled 1.940874e-05 3.221430e-06 -3.527636e-07
## DepComp_scaled:SC9_scaled
## (Intercept) 1.940874e-05
## DepComp_scaled 3.221430e-06
## SC9_scaled -3.527636e-07
## DepComp_scaled:SC9_scaled 3.156309e-04summary(interact_mod)##
## Call:
## lm(formula = PAF_scaled ~ DepComp_scaled * SC9_scaled, data = All_Variables)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5292 -0.6364 0.0136 0.6462 2.3991
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0008702 0.0159762 -0.054 0.9566
## DepComp_scaled -0.1121783 0.0169265 -6.627 3.91e-11 ***
## SC9_scaled 0.2093617 0.0159416 13.133 < 2e-16 ***
## DepComp_scaled:SC9_scaled -0.0375179 0.0177660 -2.112 0.0348 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9714 on 3710 degrees of freedom
## (1184 observations deleted due to missingness)
## Multiple R-squared: 0.05903, Adjusted R-squared: 0.05827
## F-statistic: 77.58 on 3 and 3710 DF, p-value: < 2.2e-16psych::describe(All_Variables$SC9_scaled)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 3730 0 1 0.07 0.08 1 -3.25 1.45 4.69 -0.64 0.14 0.02psych::describe(All_Variables$DepComp_scaled)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 4561 0 1 -0.16 -0.09 0.87 -2.79 7.07 9.85 1.19 3.17 0.01# Checking some model assumptions
mean(interact_mod$residuals)## [1] 1.263452e-17t.test(interact_mod$residuals)##
## One Sample t-test
##
## data: interact_mod$residuals
## t = 7.9301e-16, df = 3713, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.03123713 0.03123713
## sample estimates:
## mean of x
## 1.263452e-17plot(interact_mod)
library(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectsred = studres(interact_mod)
hist(sred, freq=FALSE,
main="Distribution of Studentized Residuals")
xfit<-seq(min(sred),max(sred),length=40)
yfit<-dnorm(xfit)
lines(xfit, yfit)
Preacher Output with School Connectedness at age 9 as the moderator
http://www.quantpsy.org/interact/mlr2.htm
TWO-WAY INTERACTION SIMPLE SLOPES OUTPUTInput
X1 = -1.58 X2 = 4.02 cv1 = -0.8 cv2 = -0.05 cv3 = 0.7 Intercept = -0.01456 X Slope = -0.17225 Z Slope = 0.2438 XZ Slope = -0.07676 df = 3700 alpha = 0.05
Asymptotic (Co)variances
var(b0) 0.0001871 var(b1) 0.00065176 var(b2) 0.00034353 var(b3) 0.00132137 cov(b2,b0) 0.0000155 cov(b3,b1) 0.00007096
Region of Significance
Z at lower bound of region = -30.775 Z at upper bound of region = -1.0859 (simple slopes are significant outside this region.)
Simple Intercepts and Slopes at Conditional Values of Z
At Z = cv1… simple intercept = -0.2096(0.0195), t=-10.7218, p=0 simple slope = -0.1108(0.0372), t=-2.9796, p=0.0029 At Z = cv2… simple intercept = -0.0267(0.0137), t=-1.9593, p=0.0502 simple slope = -0.1684(0.0255), t=-6.616, p=0 At Z = cv3… simple intercept = 0.1561(0.0194), t=8.0383, p=0 simple slope = -0.226(0.0374), t=-6.0427, p=0
Simple Intercepts and Slopes at Region Boundaries
Lower Bound…
simple intercept = -7.5175(0.5697), t=-13.1948, p=0 simple slope = 2.19(1.117), t=1.9606, p=0.05 Upper Bound…
simple intercept = -0.2793(0.0236), t=-11.8182, p=0 simple slope = -0.0889(0.0453), t=-1.9606, p=0.05
Points to Plot
Line for cv1: From {X=-1.58, Y=-0.0345} to {X=4.02, Y=-0.6552} Line for cv2: From {X=-1.58, Y=0.2393} to {X=4.02, Y=-0.7038} Line for cv3: From {X=-1.58, Y=0.5132} to {X=4.02, Y=-0.7523}
xx <- c(-1.58,4.02) # <-- change to alter plot dims
yy <- c(-1.0054,0.5132) # <-- change to alter plot dims
leg <- c(-1.58,-0.8156) # <-- change to alter legend location
x <- c(-1.58,4.02) # <-- x-coords for lines
y1 <- c(-0.0345,-0.6552)
y2 <- c(0.2393,-0.7038)
y3 <- c(0.5132,-0.7523)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='Social Deprivation',ylab='Social Deprivation')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('CVz1(1)','CVz1(2)','CVz1(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))
z1=-10 #supply lower bound for z
z2=10 #supply upper bound for z
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.17225+-0.07676*z)+(1.9606*sqrt(0.0006517551+(2*z*0.0000709585)+((z^2)*0.001321372)))
y2 <- (-0.17225+-0.07676*z)-(1.9606*sqrt(0.0006517551+(2*z*0.0000709585)+((z^2)*0.001321372)))
fy <- c(y1,y2)
fline <- (-0.17225+-0.07676*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='Moderator',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-30.775,col=4,lty=2)
abline(v=-1.0859,col=4,lty=2)
Preacher Output with Z-Scored/Centered Variables- SC9 as Moderator
TWO-WAY INTERACTION SIMPLE SLOPES OUTPUTYour Input
X1 = -2.79 X2 = 7.07 cv1 = -1 cv2 = 0 cv3 = 1 Intercept = -0.0008702 X Slope = -0.1121783 Z Slope = 0.2093617 XZ Slope = -0.0375179 df = 3710 alpha = 0.05
Asymptotic (Co)variances
var(b0) 0.00025524 var(b1) 0.0002865 var(b2) 0.00025413 var(b3) 0.00031563 cov(b2,b0) -4e-8 cov(b3,b1) 0.00000322
Region of Significance
Z at lower bound of region = -41.772 Z at upper bound of region = -1.4145 (simple slopes are significant outside this region.)
Simple Intercepts and Slopes at Conditional Values of Z
At Z = cv1… simple intercept = -0.2102(0.0226), t=-9.3142, p=0 simple slope = -0.0747(0.0244), t=-3.059, p=0.0022 At Z = cv2… simple intercept = -0.0009(0.016), t=-0.0545, p=0.9566 simple slope = -0.1122(0.0169), t=-6.6274, p=0 At Z = cv3… simple intercept = 0.2085(0.0226), t=9.2385, p=0 simple slope = -0.1497(0.0247), t=-6.0681, p=0
Simple Intercepts and Slopes at Region Boundaries
Lower Bound…
simple intercept = -8.7463(0.6661), t=-13.1305, p=0 simple slope = 1.455(0.7421), t=1.9606, p=0.05 Upper Bound…
simple intercept = -0.297(0.0276), t=-10.7468, p=0 simple slope = -0.0591(0.0301), t=-1.9606, p=0.05
Points to Plot
Line for cv1: From {X=-2.79, Y=-0.0019} to {X=7.07, Y=-0.7381} Line for cv2: From {X=-2.79, Y=0.3121} to {X=7.07, Y=-0.794} Line for cv3: From {X=-2.79, Y=0.6261} to {X=7.07, Y=-0.8499}
Who Does the Interaction NOT include
According to the Johnson-Neyman intervals, when SC9 <-1.41 stdev and when social deprivation was > 2.83 stdev, they were outside of the range of the interaction. Let’s see who that applies to. .
All_Variables %>%
filter(SC9_scaled < -1.41) %>%
summarize(n(),
mean(DepComp_scaled, na.rm = TRUE),
min(DepComp_scaled),
max(DepComp_scaled),
mean(PAF_scaled, na.rm = TRUE),
min(PAF_scaled),
max(PAF_scaled))## # A tibble: 1 x 7
## `n()` `mean(DepComp_s… `min(DepComp_sc… `max(DepComp_sc… `mean(PAF_scale…
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 341 0.100 -1.61 4.39 -0.313
## # … with 2 more variables: `min(PAF_scaled)` <dbl>, `max(PAF_scaled)` <dbl>All_Variables %>%
filter(DepComp_scaled > 2.83) %>%
summarize(n(),
mean(SC9_scaled, na.rm = TRUE),
min(SC9_scaled, na.rm = TRUE),
max(SC9_scaled, na.rm=TRUE),
mean(PAF_scaled, na.rm = TRUE),
min(PAF_scaled, na.rm = TRUE),
max(PAF_scaled, na.rm = TRUE))## # A tibble: 1 x 7
## `n()` `mean(SC9_scale… `min(SC9_scaled… `max(SC9_scaled… `mean(PAF_scale…
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 70 -0.0269 -3.07 1.28 -0.229
## # … with 2 more variables: `min(PAF_scaled, na.rm = TRUE)` <dbl>,
## # `max(PAF_scaled, na.rm = TRUE)` <dbl>All_Variables %>%
filter(DepComp_scaled > 2.83 & SC9_scaled < -1.41) %>%
summarize(n(),
mean(PAF_scaled, na.rm = TRUE),
min(PAF_scaled, na.rm = TRUE),
max(PAF_scaled, na.rm = TRUE))## # A tibble: 1 x 4
## `n()` `mean(PAF_scaled, na.rm… `min(PAF_scaled, na.rm… `max(PAF_scaled, na.rm…
## <int> <dbl> <dbl> <dbl>
## 1 3 -0.902 -1.56 0.0526Testing if the SDxSC9 Interaction Improves Fit
I calculated this based on this MPlus page: https://www.statmodel.com/chidiff.shtml
It references this paper: Satorra, A., & Bentler, P. M. (2010). Ensuring Positiveness of the Scaled Difference Chi-square Test Statistic. Psychometrika, 75(2), 243–248. https://doi.org/10.1007/s11336-009-9135-y
To Compute the Satorra-Bentler Scaled Chi-Square Difference Test:
To calculate the scaling factor:
cd=(p0 * c0 - p1*c1)/(p0 - p1) Where p=parameters and c=correction factor
p0=105
p1=107
c0=1.0297
c1=1.0298
cd=(p0*c0-p1*c1)/(p0-p1)
cd ## [1] 1.03505Then, the test statistic is as follows TRd = -2*(L0 - L1)/cd Where L=loglikelihood
L0=-65863.626
L1=-65859.958
TRd=-2*(L0-L1)/cd
TRd## [1] 7.08758This is compared to the chi-square distribution based on the degrees of freedom difference between the two models, in this case, it is 2. Therefore, our model with the interactions fits better than the main effects models.
qchisq(.95, 2)## [1] 5.991465Here is the p value for that difference test.
pchisq(TRd, df=2, lower.tail=FALSE)## [1] 0.02890357Robustness Check of the Interaction
We tested if the method for extracting the factor scores made a difference in the interpretation of the interaction. The analyses below used factor scores extracted from CFAs containing ONLY one latent variable at a time. I’m using this as the robustness check rather than the preferred method because I think the overall latent variable moderation model accounts for correlations between the latent variables and these are picked up when extracting factor scores from CFAs with all of the latent variables included in the moderation model.
All_Variables2 = All_Variables2 %>%
mutate(PAF_scaled = scale(PAF),
DepComp_scaled = scale(DepComp),
SC9_scaled = scale(SC9))
interact_mod2 = lm(PAF_scaled ~ DepComp_scaled * SC9_scaled, data = All_Variables2)
johnson_neyman(interact_mod2, pred = DepComp_scaled, modx = SC9_scaled)## JOHNSON-NEYMAN INTERVAL
##
## When SC9_scaled is OUTSIDE the interval [-35.54, -1.29], the slope of
## DepComp_scaled is p < .05.
##
## Note: The range of observed values of SC9_scaled is [-2.88, 1.18]
johnson_neyman(interact_mod2, pred = SC9_scaled, modx = DepComp_scaled)## JOHNSON-NEYMAN INTERVAL
##
## When DepComp_scaled is OUTSIDE the interval [1.23, 34.33], the slope of
## SC9_scaled is p < .05.
##
## Note: The range of observed values of DepComp_scaled is [-1.80, 7.07]
sim_slopes(interact_mod2, pred = DepComp_scaled, modx = SC9_scaled, data = All_Variables2, digits = 3)## JOHNSON-NEYMAN INTERVAL
##
## When SC9_scaled is OUTSIDE the interval [-35.540, -1.289], the slope of
## DepComp_scaled is p < .05.
##
## Note: The range of observed values of SC9_scaled is [-2.883, 1.184]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of DepComp_scaled when SC9_scaled = -0.999 (- 1 SD):
##
## Est. S.E. t val. p
## -------- ------- -------- -------
## -0.074 0.027 -2.720 0.007
##
## Slope of DepComp_scaled when SC9_scaled = -0.001 (Mean):
##
## Est. S.E. t val. p
## -------- ------- -------- -------
## -0.116 0.019 -5.986 0.000
##
## Slope of DepComp_scaled when SC9_scaled = 0.998 (+ 1 SD):
##
## Est. S.E. t val. p
## -------- ------- -------- -------
## -0.158 0.028 -5.603 0.000interact_plot(interact_mod2, pred = DepComp_scaled, modx = SC9_scaled, interval = TRUE, modx.labels = c("Low","Mean","High"), x.label = "Social Deprivation", y.label = "Positive Adolescent Function", legend.main = "School Connectedness", data = All_Variables2) + theme_classic() + theme(text=element_text(size = 16, family="serif"))
Here we check some model assumptions quickly and get the variance/covariance information necessary for the Preacher interaction tool.
# Testing assumptions about model residuals
plot(interact_mod2)
t.test(interact_mod2$residuals)##
## One Sample t-test
##
## data: interact_mod2$residuals
## t = 2.6363e-15, df = 3039, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.03490161 0.03490161
## sample estimates:
## mean of x
## 4.692634e-17mean(interact_mod2$residuals)## [1] 4.692634e-17# Getting coefficients, variances, and covariances for the Preacher interaction tool.
vcov(interact_mod2)## (Intercept) DepComp_scaled SC9_scaled
## (Intercept) 3.183108e-04 8.309214e-08 -1.264274e-08
## DepComp_scaled 8.309214e-08 3.747113e-04 2.024488e-05
## SC9_scaled -1.264274e-08 2.024488e-05 3.194435e-04
## DepComp_scaled:SC9_scaled 2.124037e-05 1.442243e-05 -3.081361e-06
## DepComp_scaled:SC9_scaled
## (Intercept) 2.124037e-05
## DepComp_scaled 1.442243e-05
## SC9_scaled -3.081361e-06
## DepComp_scaled:SC9_scaled 3.923452e-04summary(interact_mod2)##
## Call:
## lm(formula = PAF_scaled ~ DepComp_scaled * SC9_scaled, data = All_Variables2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4187 -0.6767 -0.0121 0.7054 2.1051
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.002872 0.017841 -0.161 0.8721
## DepComp_scaled -0.115903 0.019357 -5.987 2.38e-09 ***
## SC9_scaled 0.110939 0.017873 6.207 6.13e-10 ***
## DepComp_scaled:SC9_scaled -0.042074 0.019808 -2.124 0.0337 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9819 on 3036 degrees of freedom
## (1858 observations deleted due to missingness)
## Multiple R-squared: 0.02641, Adjusted R-squared: 0.02544
## F-statistic: 27.45 on 3 and 3036 DF, p-value: < 2.2e-16psych::describe(All_Variables2$DepComp)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 4561 0.01 0.57 -0.08 -0.04 0.49 -1.58 4.02 5.6 1.19 3.17 0.01psych::describe(All_Variables2$SC9)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 3332 -0.05 0.77 0.01 0.03 0.9 -2.28 0.86 3.14 -0.61 -0.17 0.01library(MASS)
sred = studres(interact_mod2)
hist(sred, freq=FALSE,
main="Distribution of Studentized Residuals")
xfit<-seq(min(sred),max(sred),length=40)
yfit<-dnorm(xfit)
lines(xfit, yfit)
Preacher Output with School Connectedness at age 9 as the moderator - Alternative Factor Score Extraction
TWO-WAY INTERACTION SIMPLE SLOPES OUTPUT
Your Input
X1 = -1.58 X2 = 4.02 cv1 = -0.82 cv2 = -0.05 cv3 = 0.72 Intercept = -0.01811 X Slope = -0.18542 Z Slope = 0.12811 XZ Slope = -0.08521 df = 3030 alpha = 0.05
Asymptotic (Co)variances
var(b0) 0.00025227 var(b1) 0.00092455 var(b2) 0.00042312 var(b3) 0.00160921 cov(b2,b0) 0.00002043 cov(b3,b1) 0.00012491
Region of Significance
Z at lower bound of region = -27.4825 Z at upper bound of region = -1.0443 (simple slopes are significant outside this region.)
Simple Intercepts and Slopes at Conditional Values of Z
At Z = cv1… simple intercept = -0.1232(0.0224), t=-5.49, p=0 simple slope = -0.1155(0.0424), t=-2.7222, p=0.0065 At Z = cv2… simple intercept = -0.0245(0.0159), t=-1.5465, p=0.1221 simple slope = -0.1812(0.0303), t=-5.9854, p=0 At Z = cv3… simple intercept = 0.0741(0.0224), t=3.3117, p=0.0009 simple slope = -0.2468(0.044), t=-5.6046, p=0
Simple Intercepts and Slopes at Region Boundaries
Lower Bound…
simple intercept = -3.5389(0.5645), t=-6.2686, p=0 simple slope = 2.1564(1.0998), t=1.9608, p=0.05 Upper Bound…
simple intercept = -0.1519(0.0259), t=-5.8637, p=0 simple slope = -0.0964(0.0492), t=-1.9609, p=0.05
Points to Plot
Line for cv1: From {X=-1.58, Y=0.0594} to {X=4.02, Y=-0.5877} Line for cv2: From {X=-1.58, Y=0.2617} to {X=4.02, Y=-0.7528} Line for cv3: From {X=-1.58, Y=0.464} to {X=4.02, Y=-0.9179}
Preacher Output with Alternative Factor Extraction and Scaled Variables - SC9 as moderator
TWO-WAY INTERACTION SIMPLE SLOPES OUTPUTYour Input
X1 = -1.58 X2 = 4.02 cv1 = -1 cv2 = 0 cv3 = 1 Intercept = -0.002005 X Slope = -0.116525 Z Slope = 0.111535 XZ Slope = -0.0423 df = 3710 alpha = 0.05
Asymptotic (Co)variances
var(b0) 0.00032174 var(b1) 0.00037875 var(b2) 0.00032289 var(b3) 0.00039657 cov(b2,b0) -1e-8 cov(b3,b1) 0.00001458
Region of Significance
Z at lower bound of region = -35.5029 Z at upper bound of region = -1.289 (simple slopes are significant outside this region.)
Simple Intercepts and Slopes at Conditional Values of Z
At Z = cv1… simple intercept = -0.1135(0.0254), t=-4.4718, p=0 simple slope = -0.0742(0.0273), t=-2.7173, p=0.0066 At Z = cv2… simple intercept = -0.002(0.0179), t=-0.1118, p=0.911 simple slope = -0.1165(0.0195), t=-5.9875, p=0 At Z = cv3… simple intercept = 0.1095(0.0254), t=4.3141, p=0 simple slope = -0.1588(0.0284), t=-5.5997, p=0
Simple Intercepts and Slopes at Region Boundaries
Lower Bound…
simple intercept = -3.9618(0.6382), t=-6.2077, p=0 simple slope = 1.3852(0.7065), t=1.9606, p=0.05 Upper Bound…
simple intercept = -0.1458(0.0293), t=-4.9759, p=0 simple slope = -0.062(0.0316), t=-1.9605, p=0.05
Points to Plot
Line for cv1: From {X=-1.58, Y=0.0037} to {X=4.02, Y=-0.4119} Line for cv2: From {X=-1.58, Y=0.1821} to {X=4.02, Y=-0.4704} Line for cv3: From {X=-1.58, Y=0.3605} to {X=4.02, Y=-0.5289}
Testing whether our results are driven by the PAF connectedness subscale
Based on the results below, it doesn’t seem like it’s driven by the connectedness subscale of the PAF scale, which is really good!
fs_scores_noconn = read_csv('FactorScores_NoConnPAF_071720.csv')## Parsed with column specification:
## cols(
## ff_id = col_character(),
## SC15 = col_double(),
## SC15_SE = col_double(),
## SC9 = col_double(),
## SC9_SE = col_double(),
## Internalizing = col_double(),
## Int_SE = col_double(),
## PAF_NoConn = col_double(),
## PAF_NoConn_SE = col_double()
## )fs_scores_noconn$ff_id = as.integer(fs_scores_noconn$ff_id)## Warning: NAs introduced by coercionAll_Variables_NoConn = read_csv('/Volumes/csmonk/csmonk-lab/Data/Fragile_Families_Longitudinal/Leigh_FFCWS/VE_SD_IntExtPAF_SchoolConn_090420.csv')## Parsed with column specification:
## cols(
## .default = col_double()
## )## See spec(...) for full column specifications.All_Variables_NoConn = All_Variables_NoConn %>%
left_join(fs_scores_noconn, by="ff_id")
is.na(All_Variables_NoConn) <- All_Variables_NoConn == 99All_Variables_NoConn = All_Variables_NoConn %>%
mutate(DepSC9 = SC9*DepComp,
DepSC15 = SC15*DepComp,
ThrSC9 = SC9*ThreatComp,
ThrSC15 = SC15*ThreatComp)summary(lm(PAF_NoConn ~ DepComp * SC9 + ThreatComp + SC15 + ThrSC9, data = All_Variables_NoConn))##
## Call:
## lm(formula = PAF_NoConn ~ DepComp * SC9 + ThreatComp + SC15 +
## ThrSC9, data = All_Variables_NoConn)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.64105 -0.41847 -0.00252 0.41836 2.27416
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0006137 0.0107682 -0.057 0.954555
## DepComp -0.0746840 0.0221917 -3.365 0.000772 ***
## SC9 0.0183789 0.0152627 1.204 0.228601
## ThreatComp 0.0693966 0.0218388 3.178 0.001497 **
## SC15 0.6973958 0.0145728 47.856 < 2e-16 ***
## ThrSC9 0.0106692 0.0303759 0.351 0.725431
## DepComp:SC9 -0.0891083 0.0314557 -2.833 0.004639 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6524 on 3704 degrees of freedom
## (1187 observations deleted due to missingness)
## Multiple R-squared: 0.4152, Adjusted R-squared: 0.4143
## F-statistic: 438.3 on 6 and 3704 DF, p-value: < 2.2e-16Digging into signifcant covariate effects
AA_Data = All_Variables %>%
filter(Race_AA == 1)
summary(lm(PAF ~ DepComp + ThreatComp + SC9 + SC15 + DepSC9 + DepSC15, data = AA_Data))##
## Call:
## lm(formula = PAF ~ DepComp + ThreatComp + SC9 + SC15 + DepSC9 +
## DepSC15, data = AA_Data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.85709 -0.45862 -0.01755 0.48511 2.07958
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.12281 0.01776 6.913 6.84e-12 ***
## DepComp -0.05139 0.03498 -1.469 0.1420
## ThreatComp 0.01097 0.03377 0.325 0.7454
## SC9 0.03622 0.02385 1.519 0.1290
## SC15 0.67216 0.02241 29.991 < 2e-16 ***
## DepSC9 -0.10203 0.04818 -2.117 0.0344 *
## DepSC15 -0.06527 0.04186 -1.559 0.1191
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6828 on 1585 degrees of freedom
## (9 observations deleted due to missingness)
## Multiple R-squared: 0.3972, Adjusted R-squared: 0.3949
## F-statistic: 174.1 on 6 and 1585 DF, p-value: < 2.2e-16C_Data = All_Variables %>%
filter(Race_C == 1)
summary(lm(PAF ~ DepComp + ThreatComp + SC9 + SC15 + DepSC9 + DepSC15, data = C_Data))##
## Call:
## lm(formula = PAF ~ DepComp + ThreatComp + SC9 + SC15 + DepSC9 +
## DepSC15, data = C_Data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.68249 -0.38639 0.02214 0.39071 2.47915
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.20773 0.02889 -7.191 2e-12 ***
## DepComp -0.15745 0.06025 -2.613 0.0092 **
## ThreatComp 0.06169 0.07835 0.787 0.4313
## SC9 0.05331 0.04040 1.320 0.1875
## SC15 0.81642 0.03600 22.675 <2e-16 ***
## DepSC9 0.05603 0.08530 0.657 0.5115
## DepSC15 -0.04518 0.07720 -0.585 0.5586
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5854 on 580 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.5747, Adjusted R-squared: 0.5703
## F-statistic: 130.6 on 6 and 580 DF, p-value: < 2.2e-16L_Data = All_Variables %>%
filter(Race_L == 1)
summary(lm(PAF ~ DepComp + ThreatComp + SC9 + SC15 + DepSC9 + DepSC15, data = L_Data))##
## Call:
## lm(formula = PAF ~ DepComp + ThreatComp + SC9 + SC15 + DepSC9 +
## DepSC15, data = L_Data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.13762 -0.43328 -0.00415 0.43563 2.23808
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.04922 0.02288 -2.152 0.031729 *
## DepComp -0.15048 0.04555 -3.304 0.000997 ***
## ThreatComp 0.03331 0.04226 0.788 0.430863
## SC9 0.02192 0.03280 0.668 0.504171
## SC15 0.78431 0.02931 26.763 < 2e-16 ***
## DepSC9 -0.03992 0.06896 -0.579 0.562781
## DepSC15 -0.07703 0.05172 -1.489 0.136790
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6282 on 801 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.5129, Adjusted R-squared: 0.5092
## F-statistic: 140.5 on 6 and 801 DF, p-value: < 2.2e-16All_Variables %>%
drop_na(Race_Cat) %>%
ggplot(aes(factor(Race_Cat), ThreatComp)) + geom_violin() + ylab('Violence Exposure') + xlab(' ')## Warning: Removed 19 rows containing non-finite values (stat_ydensity).
All_Variables %>%
drop_na(Race_Cat) %>%
ggplot(aes(factor(Race_Cat), DepComp)) + geom_violin() + ylab('Social Deprivation') + xlab(' ')## Warning: Removed 16 rows containing non-finite values (stat_ydensity).
summary(aov(ThreatComp~Race_Cat, data = All_Variables))## Df Sum Sq Mean Sq F value Pr(>F)
## Race_Cat 3 51.0 17.00 62.9 <2e-16 ***
## Residuals 3242 876.5 0.27
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 1652 observations deleted due to missingnesssummary(aov(DepComp~Race_Cat, data = All_Variables))## Df Sum Sq Mean Sq F value Pr(>F)
## Race_Cat 3 26.7 8.908 32.12 <2e-16 ***
## Residuals 3245 900.0 0.277
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 1649 observations deleted due to missingnessAll_Variables %>%
group_by(Race_Cat) %>%
summarize(VE_Mean = mean(ThreatComp, na.rm=TRUE), SD_Mean = mean(DepComp, na.rm = TRUE))## `summarise()` ungrouping output (override with `.groups` argument)## # A tibble: 5 x 3
## Race_Cat VE_Mean SD_Mean
## <fct> <dbl> <dbl>
## 1 African American 0.124 0.0451
## 2 Caucasian -0.192 -0.189
## 3 Latinx -0.0782 0.0549
## 4 Other 0.0235 0.00572
## 5 <NA> -0.00396 0.0132Who is the Primary Caregiver?
Age 3: IH3int5 01 Bio Mother 02 Bio Father 03 Mat Grandmother 04 Pat Grandmother 05 Other Relative 06 Other Non-Relative 07 Mat Grandmother 08 Mat Grandfather
Age5: IH5int5 01 Bio Mother 02 Bio Father 03 Mat Grandmother 04 Pat Grandmother 05 Other Relative 06 Other Non-Relative 07 Mat Grandmother 08 Mat Grandfather
Age 9: pcg5idstat 61 Bio Mother
63 Other
62 Bio Father
PC_Info = read_csv('primary_caregiver_relationships.csv')## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_double(),
## ff_id = col_double(),
## IH3int5 = col_double(),
## IH5int5 = col_double(),
## pcg5idstat = col_double()
## )PC_Info %>%
group_by(IH3int5) %>%
summarize(n())## `summarise()` ungrouping output (override with `.groups` argument)## # A tibble: 7 x 2
## IH3int5 `n()`
## <dbl> <int>
## 1 1 3191
## 2 2 15
## 3 3 10
## 4 4 3
## 5 5 4
## 6 6 5
## 7 NA 1670PC_Info %>%
group_by(IH5int5) %>%
summarize(n())## `summarise()` ungrouping output (override with `.groups` argument)## # A tibble: 7 x 2
## IH5int5 `n()`
## <dbl> <int>
## 1 1 2914
## 2 2 27
## 3 3 19
## 4 4 7
## 5 5 10
## 6 6 1
## 7 NA 1920PC_Info %>%
group_by(pcg5idstat) %>%
summarize(n())## `summarise()` ungrouping output (override with `.groups` argument)## # A tibble: 4 x 2
## pcg5idstat `n()`
## <dbl> <int>
## 1 61 3828
## 2 62 180
## 3 63 194
## 4 NA 696Correlation between Age 9 and Age 15 CBCL scores
cor.test(All_Variables$Intern_CBCL15, All_Variables$InternCBCL9)##
## Pearson's product-moment correlation
##
## data: All_Variables$Intern_CBCL15 and All_Variables$InternCBCL9
## t = 17.068, df = 2791, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2734574 0.3406337
## sample estimates:
## cor
## 0.3074286cor.test(All_Variables$Extern_CBCL15, All_Variables$ExternCBCL9)##
## Pearson's product-moment correlation
##
## data: All_Variables$Extern_CBCL15 and All_Variables$ExternCBCL9
## t = 28.541, df = 2850, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.4424257 0.4995323
## sample estimates:
## cor
## 0.4714731Extra Analyses
Regressions with externalizing outside of Mplus
Age 9
summary(lm(Ext9 ~ DepComp + ThreatComp, data = All_Variables))
summary(lm(Ext9 ~ DepComp + ThreatComp + SC9, data = All_Variables))
summary(lm(Ext9 ~ DepComp + ThreatComp + SC9 + DepSC9 + ThrSC9, data = All_Variables))Age 15
summary(lm(Ext15 ~ DepComp + ThreatComp, data = All_Variables))
summary(lm(Ext15 ~ DepComp + ThreatComp + SC9 + SC15, data = All_Variables))
summary(lm(Ext15 ~ DepComp + ThreatComp + SC9 + SC15 + DepSC9 + DepSC15 + ThrSC9 + ThrSC15, data = All_Variables))Regressions with internalizing outside of Mplus
Age 9
summary(lm(Int9 ~ DepComp + ThreatComp, data = All_Variables))
summary(lm(Int9 ~ DepComp + ThreatComp + SC9, data = All_Variables))
summary(lm(Int9 ~ DepComp + ThreatComp + SC9 + DepSC9 + ThrSC9, data = All_Variables))Age 15
summary(lm(Int ~ DepComp + ThreatComp, data = All_Variables))
summary(lm(Int ~ DepComp + ThreatComp + SC9 + SC15, data = All_Variables))
summary(lm(Int ~ DepComp + ThreatComp + SC9 + SC15 + DepSC9 + DepSC15 + ThrSC9 + ThrSC15, data = All_Variables))Regressions with PAF outside of Mplus
summary(lm(PAF ~ DepComp + ThreatComp + SC9 + SC15 + DepSC9, data = All_Variables))
summary(lm(PAF ~ DepComp + ThreatComp + SC9 + SC15 + DepSC9 + DepSC15, data = All_Variables))
summary(lm(PAF ~ DepComp + ThreatComp + SC9 + SC15 + DepSC9 + DepSC15 + cm1bsex, data = All_Variables))
summary(lm(PAF ~ DepComp + ThreatComp + SC9 + SC15 + DepSC9 + DepSC15 + cm1bsex + povco_avg + Race_AA + Race_C + Race_L, data = All_Variables))Possible Interaction with SD and SC9 for Ext9
All_Variables = All_Variables %>%
mutate(Ext9_scaled = scale(Ext9))
DepSC9 = lm(Ext9_scaled ~ DepComp_scaled * SC9_scaled + ThreatComp + ThrSC9, data = All_Variables)
summary(DepSC9)##
## Call:
## lm(formula = Ext9_scaled ~ DepComp_scaled * SC9_scaled + ThreatComp +
## ThrSC9, data = All_Variables)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0309 -0.6895 -0.0416 0.6023 5.2211
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.007491 0.016116 -0.465 0.6421
## DepComp_scaled 0.119760 0.019141 6.257 4.44e-10 ***
## SC9_scaled -0.110588 0.015366 -7.197 7.59e-13 ***
## ThreatComp 0.541467 0.033670 16.082 < 2e-16 ***
## ThrSC9 0.029420 0.043641 0.674 0.5003
## DepComp_scaled:SC9_scaled 0.045613 0.018774 2.430 0.0152 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9243 on 3305 degrees of freedom
## (1587 observations deleted due to missingness)
## Multiple R-squared: 0.1421, Adjusted R-squared: 0.1408
## F-statistic: 109.5 on 5 and 3305 DF, p-value: < 2.2e-16interact_plot(DepSC9, pred = DepComp_scaled, modx = SC9_scaled, interval = TRUE, x.label = "Social Deprivation", y.label = "Externalizing Sx Age 9", legend.main = "School Connectedness", data = All_Variables) + theme_classic() 
interact_plot(DepSC9, pred = SC9_scaled, modx = DepComp_scaled, interval = TRUE, x.label = "School Connectedness", y.label = "Externalizing Sx Age 9", legend.main = "Social Deprivation", data = All_Variables) + theme_classic() 
johnson_neyman(DepSC9, pred = DepComp_scaled, modx = SC9_scaled)## JOHNSON-NEYMAN INTERVAL
##
## When SC9_scaled is OUTSIDE the interval [-13.64, -1.31], the slope of
## DepComp_scaled is p < .05.
##
## Note: The range of observed values of SC9_scaled is [-3.25, 1.45]
johnson_neyman(DepSC9, pred = SC9_scaled, modx = DepComp_scaled)## JOHNSON-NEYMAN INTERVAL
##
## When DepComp_scaled is OUTSIDE the interval [1.22, 12.75], the slope of
## SC9_scaled is p < .05.
##
## Note: The range of observed values of DepComp_scaled is [-1.80, 7.07]
sim_slopes(DepSC9, pred = DepComp_scaled, modx = SC9_scaled, data = All_Variables, digits = 3)## JOHNSON-NEYMAN INTERVAL
##
## When SC9_scaled is OUTSIDE the interval [-13.644, -1.306], the slope of
## DepComp_scaled is p < .05.
##
## Note: The range of observed values of SC9_scaled is [-3.246, 1.445]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of DepComp_scaled when SC9_scaled = -1.055 (- 1 SD):
##
## Est. S.E. t val. p
## ------- ------- -------- -------
## 0.072 0.027 2.638 0.008
##
## Slope of DepComp_scaled when SC9_scaled = -0.005 (Mean):
##
## Est. S.E. t val. p
## ------- ------- -------- -------
## 0.120 0.019 6.246 0.000
##
## Slope of DepComp_scaled when SC9_scaled = 1.045 (+ 1 SD):
##
## Est. S.E. t val. p
## ------- ------- -------- -------
## 0.167 0.028 6.024 0.000Based on the simple slopes/regions of significance, it seems like this interaction covers a similar story as the one predicting PAF, but in a much smaller range of values. It seems like school connectedness is more protective for those with lower social deprivation, but the regions wehre that is the case are much smaller and I think it will not add to the main effects that we have which show the same thing.